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Abstract.  We  begin  this  paper  by  identifying  a  class  of  stochastic  mixed-integer  programs  that  have  col¬ 
umn-oriented  formulations  suitable  for  solution  by  a  branch-and-price  algorithm  (B&P).  We  then  survey  a 
number  of  examples,  and  use  a  stochastic  facility-location  problem  (SFLP)  for  a  detailed  demonstration  of  the 
relevant  modeling  and  solution  techniques.  Computational  results  with  a  scenario  representation  of  uncertain 
costs,  demands  and  capacities  show  that  B&P  can  be  orders  of  magnitude  faster  than  solving  the  standard 
formulation  by  branch  and  bound.  We  also  demonstrate  how  B&P  can  solve  SFLP  exactly  -  as  exactly  as 
a  deterministic  mixed-integer  program  -  when  demands  and  other  parameters  can  be  represented  as  certain 
types  of  independent,  random  variables,  e.g.,  independent,  normal  random  variables  with  integer  means  and 
variances. 
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1.  Introduction 

This  paper  describes  a  class  of  two-stage  stochastic  mixed-integer  programs  (SMIPs) 
whose  instances  are  amenable  to  column-oriented  formulations,  and  then  shows  how  to 
solve  such  formulations  with  a  branch-and-price  algorithm  (B&P).  The  phrase  “branch 
and  price”  was  coined  by  Savelsbergh  [49],  but  the  technique  was  first  proposed  by  John¬ 
son  [33]  and  implemented  by  Desrochers  and  Solomon  [25]  and  Desrochers  et  al.  [24]. 
Branch  and  price  combines  dynamic  column  generation  -  this  is  known  widely  through 
the  “cutting-stock  problem”  of  Gilmore  and  Gomory  [30]  -  with  standard  branch  and 
bound. 

Stochastic  programmers  have  only  just  begun  to  see  that  B&P  applies  to  their  prob¬ 
lems,  and  we  find  only  two  papers  on  the  topic:  Damodaran  and  Wilhelm  [21]  and  Lulli 
and  Sen  [42],  (However,  Shiina  and  Birge  [51]  and  Singh  et  al.  [53]  use  column-gen¬ 
eration  without  branch  and  bound  to  solve  SMIPs.)  Those  papers  investigate  specific 
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applications  of  B&P  to  stochastic  programming.  In  contrast,  we  describe  a  complete 
class  of  problems  to  which  B&P  applies;  in  similarity,  we  show  impressive  computa¬ 
tional  results. 

After  defining  the  special  class  of  SMIPs,  we  provide  examples  to  show  how  sto¬ 
chastic  versions  of  several  well-known  deterministic  models  fit  into  this  framework: 
the  elastic  generalized-assignment  problem  (Brown  and  Graves  [12]),  crew-scheduling 
(Vance  et  al.  [56],  Day  and  Ryan  [23]),  vehicle-routing  problems  (e.g.,  Desrosiers  et  al. 
[26]),  and  the  origin-destination  integer  multicommodity  flow  problem  (Barnhart  et  al. 
[6]).  One  additional  problem,  a  stochastic  facility-location  problem  (SFLP),  guides  our 
detailed  exploration  of  the  B&P  solution  approach. 

We  initially  model  and  solve  a  version  of  SFLP  with  uncertain  demands,  costs  and 
capacities,  all  represented  through  scenarios.  Such  representations  appear  frequently 
in  the  stochastic-programming  literature  (e.g.,  Butler  and  Dyer  [15],  Chen  et  al.  [18], 
Ahmed  and  Sahinidis  [1],  Lulli  and  Sen  [42]),  with  the  primary  advantage  being  to  allow 
arbitrary  dependence  among  uncertain  parameters.  However,  another  common  formu¬ 
lation  approach  defines  individual  probability  distributions  for  the  stochastic  program’s 
parameters,  which  are  typically  assumed  to  be  independent  (e.g.,  Bertsimas  [  10],  Zhou 
and  Liu  [61]).  We  will  show  how  B&P  can  solve  SFLP  in  this  situation,  too.  Further¬ 
more,  we  will  solve  the  problem  exactly,  that  is,  with  the  same  certitude  that  prevails 
in  deterministic  mixed-integer  programs.  This  contrasts  with  alternative  solution  proce¬ 
dures  that  provide  only  probabilistic  or  asymptotic  guarantees  of  solution  quality  (e.g., 
Carpe  and  Tind  [17],  Sen  and  Higle  [50],  Ahmed  and  Sahinidis  [1]). 

Solution  methods  for  SMIPs  with  “scenario  uncertainty’’  typically  employ  Benders 
decomposition  (Benders  [9]),  or  extensions  thereof,  to  the  original  model.  Examples 
include  the  integer  L-shaped  decomposition  from  Laporte  and  Louveaux  [36],  and  the 
methods  developed  by  Carpe  and  Tind  [17]  and  by  Sen  and  Higle  [50].  Unfortunately, 
all  of  these  decompositions  use  a  master  problem  whose  linear-programming  relaxation 
is  no  stronger  than  the  linear-programming  relaxation  of  the  original  model  (measured 
over  the  master-problem  variables  which  the  two  formulations  have  in  common).  Con¬ 
sequently,  these  decompositions  will  suffer  if  the  original  model  formulation  has  a  poor 
continuous  relaxation.  In  contrast,  B&P  solves  a  column-oriented  reformulation  of  a 
model,  also  by  a  form  of  decomposition,  but  that  reformulation  will  normally  have  a 
tighter  relaxation  than  the  original  model  (Barnhart  et  al.  [7]). 

The  next  section  describes  our  special  class  of  SMIPs  and  shows  how  to  convert 
their  standard  formulations  into  column-oriented  ones.  Several  models  from  the  litera¬ 
ture  provide  concrete  examples.  Section  3  presents  the  SFLP  with  scenario  uncertainty, 
describes  how  to  solve  instances  with  B&P,  and  provides  computational  results.  Section 
4  presents  a  version  of  SFLP  in  which  random  parameters  take  on  continuous  distribu¬ 
tions,  describes  how  to  solve  instances  with  B&P,  and  provides  computational  results. 
Section  5  presents  conclusions. 


2.  General  Methodology 

We  will  show  that  a  variety  of  SMIPs  of  the  following  form  can  arise  in  scheduling, 
routing  and  other  applications: 
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Formulation  (SMIPO): 


min  x,  +  (1) 

iel 

s.t.y^  AjXj  =  b  (2) 

iel 

x,  e  Xt  Vi  e  /  (3) 

where,  for  all  iel, 

hi(xi,  f(  )  =  min  f(-y,-  (4) 

yt 

s.t.  Ay,  >  BjXj  +  d,  (5) 

y,  e  Yi,  (6) 


and  where  =  vec (A, ,  Dj ,  d, ,  ?,  ).  The  sets  X,  require  all  x,  to  be  bounded  and  integral. 
The  sets  T,  will  normally  require  non-negativity  of  the  y,-,  but  we  place  no  particular 
restrictions  on  these  sets;  the  variables  may  be  continuous  and/or  integer.  The  objective- 
function  term  ^  A.  [/z,-  (x,- ,  |,  )]  is  called  the  recourse  function  (Walkup  and  Wets  [57]). 
iel  ' 

We  further  assume  that  the  model  exhibits  relatively  complete  recourse  (Rockafellar 
and  Wets  [47]),  which  implies  that,  for  any  x;  e  X,,  an  optimal  solution  y,-,  satisfying 
constraints  (5)  and  (6),  can  always  be  found.  We  note  that  D ,  is  an  identity  matrix  in 
some  of  our  examples,  implying  the  property  of  simple  recourse  (Beale  [8],  Wets  [59]). 
However,  this  is  not  an  inherent  requirement  of  this  class  of  problems.  The  key  feature 
of  SMIPO  is  that  the  recourse  function  decomposes  by  “subproblem”  i. 

Because  all  x,  are  integral  and  bounded,  the  principles  of  Dantzig- Wolfe  decompo¬ 
sition  apply  (Dantzig  and  Wolfe  [22]),  as  extended  to  integer  programs  by  Appelgren 
[4].  (See  Wolsey  [60],  Section  11.2,  for  a  comprehensive  discussion.)  To  describe  this 
decomposition,  let  x*  e  X;,k  e  Kj ,  denote  the  enumerated  vectors  satisfying  con¬ 
straints  (3)  for  each  /.  Because  of  relatively  complete  recourse,  E~.  [/;,  (x^\  f, )]  is  well 
defined  for  all  such  xi .  Then,  because  of  the  special  structure,  we  can  embed  the  decom¬ 
posed  recourse  function  into  the  column  costs  of  a  column-oriented  formulation  for 
SMIPO: 

Column-Oriented,  Mixed-Integer,  Two-Stage  Stochastic  Program  (CSMIPO) 
Indices: 

i  e  I  subproblems 

k  e  Kj  indices  for  (integer)  vectors  x^  e  X ; 
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Data: 

xf  the  Alh  vector  xf  e  X, 

Cj  first-stage  costs  for  subproblem  i 

Decision  Variables: 

A.f  1  if  xf  e  Xi  is  selected  as  subproblem  ;’s  solution  in  the  overall  solution,  and  0 
otherwise 

Formulation  (CSMIPO): 

min  (c'*f  +  Eh  [*(*f  •  fi)])  (7) 

iel  keKj 


8-tEEMi/=b  <8> 

iel  keKj 

^lf=l  V/e/  (9) 

keKj 

Xk  e  {0,  1}  V  i,  k  (10) 

Constraints  (9)  are  often  referred  to  as  “convexity  constraints.”  CSMIPO  can  be  applied 
when  first-stage  variables  are  general  integers,  but  binary  variables  are  typical  so  we 
assume  this  restriction  hereafter  for  simplicity. 

Two  examples  of  CSMIPO  to  be  discussed  shortly  simplify  to  the  form  of  a  set-par¬ 
titioning  model,  so  we  describe  this  special  case  first.  If  constraints  (8)  are  indexed  by 
j  e  /,  b  is  a  vector  of  Is,  A;  is  the  identity  matrix  for  all  i,  and  all  vectors  xf  are  binary, 
0-1,  then  CSMIPO  may  be  written  as: 


(11) 


s-1-  EE^‘  =  1  vJeJ 

(12) 

iel  keKj 

M 

II 

< 

m 

(13) 

keKj 


x)  e  {0,  1}  V;  el,  k  e  Kt 
where  xfj  =  (A,xf) j  and  cf  =  c ,xf  +  £|.  [A (if ,  f,)]. 


(14) 
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Naturally,  the  cardinalities  of  the  index  sets  K,  in  CSMIPO  or  CSMIP1  may  be 
enormous,  and  it  will  usually  be  necessary  to  solve  these  models  without  explicitly 
enumerating  the  x^ .  Before  discussing  such  issues,  however,  we  will  provide  concrete 
examples  of  how  this  reformulation  technique  applies  to  some  stochastic  versions  of 
well-known,  deterministic  optimization  problems.  We  supply  only  short  descriptions  of 
the  problems,  and  ask  the  reader  to  check  the  references  for  more  details.  For  simplicity, 
we  hereafter  drop  the  subscript  on  the  expectation  operator,  because  it  should  be  clear 
from  context. 


2.1.  Elastic  generalized-assignment  problem 
(Brown  and  Graves  [12],  Appleget  and  Wood  [3]) 

The  objective  of  the  (deterministic)  elastic  generalized-assignment  problem  (EGAP)  is 
to  minimize  the  cost  of  assigning  capacity-consuming  tasks  j  e  J  to  capacitated  agents 
i  e  /,  so  that  (a)  each  task  is  assigned  to  exactly  one  agent,  and  (b)  the  total  capacity 
assigned  to  agent  i  does  not  exceed  the  agent’s  (potentially  uncertain)  capacity  hj ,  unless 
an  appropriate  per-unit  penalty  f)  is  paid.  If  assigning  task  j  to  agent  i  consumes  a  random 
amount  of  agent  i’s  capacity,  denoted  bjj,  and  if  we  denote  the  direct  cost  of  such  an 
assignment  as  Cjj  (which  could  represent  an  expected  value),  a  stochastic  version  of  the 
EGAP  (SEGAP)  can  be  defined  (Spoerl  and  Wood  [54]): 

SEGAP 

X>»,+  £[ftj(x^(bj,  «,•))])  (15) 

i€l  \jsJ  ) 

s.t.  J2x‘j  ~  '  ^ j  £  J  (16) 

iel 

Xij  e  [0,  1}  Vi  e  /,  j  e  J  (17) 

where 

hi  (x; ,  (b,-,  hi))  =  min yi  (18) 

yi 

s-t.  Vi  >  ^2  bj  j  xij  -  hi  (19) 

j£j 

yt  >  0.  (20) 

Here,  xtj  equals  1  if  task  j  is  assigned  to  agent  i  and  is  0  otherwise,  and  y;  represents 
capacity  violation  for  agent  i.  Therefore,  li [/? ,( x, ,  (b;,  m,))]  represents  the  expected 
capacity- violation  penalty  for  agent  i. 

The  conversion  of  SEGAP  to  a  column-oriented  formulation  is  straightforward.  (Sa- 
velsbergh  [49]  creates  the  analogous  formulation  for  a  deterministic,  inelastic  GAP.) 
Each  variable  represents  a  potential  joint  assignment  of  tasks  to  a  particular  agent,  i.e., 
a  collection  of  tasks  that  an  agent  might  be  required  to  perform. 
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Column-Oriented  Formulation  for  SEGAP  (CSEGAP) 


Indices: 

i  e  I  agents 
j  e  J  tasks 

A  e  Kj  joint  assignments  of  tasks  to  agent  i 

Data: 

xfj  1  if  task  j  is  assigned  to  agent  i  in  the  Ath  joint  assignment  for  agent  i,  and  0 
otherwise 

xf  . . .  ,  the  Ath  joint- assignment  vector  for  agent  i 

Xj  the  set  of  all  possible  joint  assignments  for  agent  i  (the  index  set  Kj  can  now 
be  defined  as  the  minimal  set  such  that  U^gj^xf  =  Xj ) 
cf  expected  cost  ofthe  Ath  joint  assignment  for  agent  t'(cf  =  CjXl>+E[hj(xll ,  (b,-,  it,-))] 

for  all  i  e  I  and  k  e  Kj ) 

Decision  Variables: 

A*  1  if  the  Ath  joint  assignment  is  selected  for  agent  i,  and  0  otherwise 

Formulation  (CSEGAP):  Same  as  CSMIP1,  Equations  ( 1 1 )— ( 1 4) . 

Constraints  (12)  guarantee  that  each  task  j  is  assigned  to  exactly  one  agent,  and  con¬ 
straints  (13)  ensure  each  agent  i  receives  exactly  one  joint  assignment  of  tasks.  Note  that 
this  simple  model  allows  any  group  of  tasks  to  be  assigned  to  any  agent,  so  Xj  consists 
of  all  binary  vectors  of  length  |/|,  implying  that  K,  =  2^  for  all  i.  This  model  can, 
indeed,  possess  a  large  number  of  columns. 


2.2.  Routing  and  scheduling  with  time  windows 
(Desrosiers  et  al.  [26],  Ribeiro  and  Soumis  [46]) 

The  vehicle  routing  problem  with  time  windows  (VRPTW)  is  one  important  exemplar 
from  our  special  problem  class.  VRPTW  describes  a  fleet  of  vehicles  that  must  deliver 
a  set  of  customer  orders,  with  each  customer  being  represented  by  a  node  in  a  network. 
Vehicles  have  limited  capacities,  and  customers  specify  windows  of  time  during  which 
deliveries  should  be  made.  The  model  allocates  orders  to  vehicles  so  that  vehicle  capac¬ 
ities  are  respected,  and  identifies  a  route  for  each  vehicle  that  delivers  each  order  during 
that  order’s  customer-specified  time  window.  The  column-oriented  formulation  for  this 
problem  fits  the  form  of  CSMIP1,  Equations  (1 1)— (14),  with  indices,  parameters  and 
variables  appropriately  defined.  The  parameters  x*,  A  e  Kj,  now  represent  potential 
routes  and  sets  of  deliveries  to  customers  for  vehicle  ;,  each  of  which  covers  a  subset  of 
the  customer  set  J.  A  probabilistic  recourse  function  could  include  expected  penalties  for 
violating  time  windows  and  expected  penalties  for  exceeding  maximum  route  duration. 
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2.3.  Crew  scheduling 

(Vance  et  al.  [56],  Day  and  Ryan  [23]) 

Following  the  description  in  [56],  an  airline  crew-scheduler  wishes  to  minimize  the  cost 
of  assigning  flight  crews  to  a  fixed  schedule  of  flights.  Crew  pairings  define  feasible  trip 
itineraries  that  can  be  assigned  to  some  crew.  Each  pairing  consists  of  a  sequence  of 
flights  that  starts  and  ends  at  a  home  base,  respects  limits  on  work  hours,  allows  times 
for  rest  breaks,  and  satisfies  numerous  other  requirements.  Set-partitioning  models  are 
the  norm  for  this  type  of  problem  (e.g.,  [56],  [23]),  and  these  fit  a  simplified  form  of 
CSMIPO  in  which  the  set  /  is  a  singleton,  b  is  a  vector  of  Is  corresponding  to  flights  that 
must  be  covered  by  crews,  and  the  convexity  constraints  (9)  are  eliminated.  The  binary 
columns  A,  xf  represent  pairings. 

However,  “home-base  constraints”  may  need  to  be  enforced  (e.g.,  Butchers  et  al. 
[14])  and  these  simply  modify  constraints  (9)  to 

^  11  i  V/  e  7’  (21) 

ksKi 

where  I  denotes  the  set  of  home  bases,  uj  denotes  the  number  of  crews  available  at 
i  e  /,  and  the  index  sets  K,  now  represents  potential  pairings  for  crews  based  at  i.  For 
the  recourse  function,  we  suggest  a  probabilistic  variant  on  the  function  that  Ehrgott 
and  Ryan  [28]  use  to  penalize  schedules  that  do  not  allow  adequate  time  for  crews  to 
switch  aircraft.  Their  function  is  based  on  averaged  historical  information,  but  could 
be  modified  to  represent  a  penalty  function  integrated  over  empirical  or  fitted  delay 
distributions. 


2.4.  Origin-destination  integer  multi-commodity  flow  problem 
(Barnhart  et  al.  [6]) 

The  origin-destination  integer  multi-commodity  flow  problem  restricts  the  standard,  lin¬ 
ear  multi-commodity  flow  problem  (Ahuja  et  al.  [2],  chapter  17)  by  requiring  each  of  a 
set  of  commodities  to  be  shipped  from  its  origin  to  its  destination  along  a  single  path.  This 
problem’s  formulation  resembles  the  well-known  path-oriented  (i.e.,  column-oriented) 
formulation  of  the  linear  multi-commodity  flow  problem  (Ford  and  Fulkerson  [29]),  but 
with  binary  variables.  In  particular,  a*  =  1  if  commodity  i  follows  path  k  £  Kj,  and 
)J-  —  0  otherwise.  The  main  constraints  of  this  problem  require  that  the  sum  of  all  com¬ 
modities  flowing  across  each  arc  respect  that  arc’s  capacity.  Constraints  (8),  converted  to 
inequalities,  handle  these  requirements  if  we  (a)  let  bj  represent  the  capacity  of  arc  /,  and 
(b)  define  ctijxf  to  be  the  amount  of  arc  j’s  capacity  consumed  if  commodity  i  is  shipped 
using  path  k.  The  convexity  constraints  (9)  guarantee  selection  of  a  single  path,  with 
appropriate  origin  and  destination,  for  each  commodity.  For  a  communications  network, 
say,  each  component  of  the  recourse  function  E  [ h ,  (x,- ,  f ,• )]  might  represent  an  expected, 
path-dependent  penalty  based  on  uncertain  link  availability  (Girard  and  Sanso  [31])  or 
uncertain  “hop  delay”  that  is  independent  of  congestion  (Papagiannaki  et  al.  [44]). 

The  following  section  investigates,  in  detail,  one  additional  problem  that  fits  the 
framework  of  SMIPO,  CSMIPO,  and  more  specifically,  CSMIP1. 
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3.  Solving  a  stochastic  facility  location  problem  by  branch  and  price 

3.1.  A  stochastic  facility  location  problem  with  sole  sourcing 


A  standard,  deterministic,  facility-location  problem  seeks  the  best  locations  at  which  to 
open,  i.e.,  build  or  lease,  capacitated  production  facilities  that  will  ship  to  established 
customers  to  meet  those  customers*  demands  for  some  product.  (We  consider  only  the 
simplest,  single -product  case  here.)  The  mathematical  model  must  find  the  best  trade¬ 
off  between  variable  and  fixed  costs  (Laporte  et  al.  [36]):  Opening  more  facilities  leads 
to  lower  shipping  (variable)  costs  because  plants  are  closer  to  customers,  on  average, 
but  opening  more  facilities  incurs  more  facility-installation  (fixed)  costs.  The  determin¬ 
istic  model  typically  assumes  that  all  customer  demands  will  be  completely  satisfied, 
and  sometimes  requires  that  each  customer  be  served  by  a  unique  facility.  This  latter 
requirement  is  known  as  sole-sourcing ,  and  the  resulting  model  is  called  the  (deter¬ 
ministic)  capacitated  facility-location  problem  with  sole-sourcing  (FLP)  (Barcelo  and 
Casanova  [5]). 

Assume  now  that  some  uncertainty  arises  in  the  nominally  deterministic  FLP:  Does 
a  manufacturer  really  know  what  future  costs,  demands  and  capacities  will  be?  Let  us, 
initially,  represent  this  uncertainty  through  a  finite,  discrete  set  of  scenarios  indexed  by  s, 
with  probability  of  occurrence  equaling  ps ,  and  with  c'  ,  u ' ,  d'  and  f  '  representing  cor¬ 
responding  shipping  costs,  facility  capacities,  customer  demands  and  “excess-demand” 
penalties,  respectively.  For  simplicity,  we  assume  that  if  the  aggregate  demand  for  a 
facility  exceeds  its  capacity  to  produce,  the  facility  pays  a  penalty  based  on  the  excess. 
Thus,  the  model  that  follows  will  be  reasonable  if,  when  a  facility  runs  short  of  sup¬ 
ply,  it  purchases  extra  product  from  an  outside  supplier  and  actually  satisfies  the  excess 
demand  by  shipping  this  product  appropriately. 


Stochastic  Facility  Location  Problem  with  Sole-Sourcing  (SFLP) 
Indices: 

i  <=  I  potential  facilities  (at  various  locations) 
j  e  J  customers 
s  e  S  scenarios 


Data  [units]: 

c'  fixed  cost  to  open  facility  i  [dollars] 

dSj  customer  j’s  demand  under  scenario  s  [tons] 

w?  facility  V s  capacity  under  scenario  s  [tons] 

f?  penalty  for  each  unit  of  excess  demand  at  facility  i  under  scenario  .y  [dollars/ton] 
ps  probability  that  scenario  s  occurs 

Cfj  per-unit  shipping  cost  from  i  to  j  under  scenario  .y  [$/ton] 

Cij  expected  cost  to  supply  all  of  customer  j’s  demand  from  facility  i,  disregarding 
any  shortfalls  in  facility  f  s  capacity  [dollars]  (c,y  =  ^ 

seS  7  7 
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Decision  Variables  [units]: 

x'  1  if  facility  i  is  opened,  and  0  otherwise 

Xij  1  if  customer  j  is  assigned  to  facility  i,  and  0  otherwise 

yj  amount  of  excess  demand  at  facility  i  under  scenario  .v  [tons] 


Formulation  (SFLP): 


X  ,x,y  ,e/ 

S.t.  -  x\  +  Xij 

iel 


+  cijxij  +  J2H  psf>ysi 

iel  jeJ  seS  iel 

<  0  V/  €  I,  j  e  J 
=  1  V./  e  J 


J2 d)xu ~ yt  - ui  Vie/’seS 

jsj 

x\  e  [0,  1}  Vi  e  I 
xij  e  [0,  1}  Vi  e  /,  j  e  J 
y-  >0  Vi  e  I,  s  e  S 


(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


This  type  of  formulation  is  known  as  the  extensive  form  of  a  stochastic  program 
(Birge  and  Louveaux  [1 1],  p.  8),  because  the  second-stage  variables  and  constraints  are 
made  explicit  for  all  scenarios. 


3.2.  A  column-oriented  formulation  for  SFLP 

Here  we  describe  CSFLP,  a  column-oriented  formulation  for  SFLP  that  fits  directly  into 
the  form  of  CSMIP1.  In  this  formulation,  a  joint  assignment  represents  any  collection 
of  customers  that  are  served  by  the  same  facility.  We  note  that  Teo  and  Shu  [55]  and  Lo- 
rena  and  Senne  [38]  have  previously  used  column  generation  for  solving  deterministic 
facility-location  problems.  Both  of  these  papers  use  master  problems  that  resemble  our 
set-partitioning  master  problem.  The  main  differences  lie  in  the  subproblems:  Theirs  are 
deterministic  while  ours  is  stochastic. 


Column-Oriented  Formulation  of  SFLP  (CSFLP) 


Indices: 


k  e  Kj  possible  joint  assignments  of  customers  to  facility  i 
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Data  [units]: 

xfj  1  if  customer  j  is  assigned  to  facility  i  in  the  kth  joint  assignment  for  that  facility, 
and  0  otherwise 

cf  total  expected  cost  of  the  kth  joint  assignment  for  facility  i  (cf  —  c\  +  ^ 

+J]  Psfi9i>  except  cf  =  0  for  the  null  assignment)  [dollars] 

5 


Decision  Variables: 

/,f  1  if  the  Ath  joint  assignment  of  customers  to  facility  i  is  selected,  and  0  otherwise 
Formulation  (CSFLP):  Same  as  CSMIP1,  Equations  (11) — (14). 

A  column-oriented  formulation  like  CSFLP  cannot  be  solved  directly  because  it  is 
impossible,  or  impractical,  to  create  the  full  set  of  columns.  Therefore,  each  set  K ,  is 
replaced  by  a  subset  to  form  a  restricted  master  problem  (RMP).  The  solution  to  the 
LP  relaxation  of  the  RMP  (LP-RMP)  then  yields  dual  variables,  which  can  be  used  to 
attempt  to  identify  one  or  more  new  columns  with  favorable  reduced  costs  through  one 
or  more  column-generation  subproblems.  In  the  case  of  CSFLP,  if  we  seed  the  RMP 
with  all  null  joint  assignments,  the  following  subproblem  arises  for  each  facility  i: 

CSUB;(7T,  Pi) 


min£  (dij-Tij 

0  *ij  +  J2  +  ci  ~  & 

(29) 

jeJ 

seS 

s.t.  J2  djx<J 

-  yf  -  u\  VseS 

(30) 

jeJ 

Xij  e  [0,  1}  Vj  e  J 

(31) 

ysi  >0  VseS, 

(32) 

where  if  /  is  the  optimal  dual  variable  associated  with  constraint  (12)  for  customer  j  in 
LP-RMP,  and  p,  is  the  optimal  dual  variable  from  LP-RMP  for  the  convexity  constraint 
(13)  associated  with  facility  i.  If  the  per-unit  shipping  cost  from  facility  i  to  customer  j 
is  a  random  quantity,  denoted  c,y,  note  that  cf/  =  E[cjjdj],  and  q7  =  E[£jj]E[dj]  if 
independence  of  Cjj  and  dj  prevails. 

If  the  solution  to  CSUB;  (jf ,  Pi )  defines  a  non-null  joint  assignment  of  customers  to 
facility  i,  z*  gives  the  reduced  cost  of  the  joint  assignment  with  respect  to  LP-RMP’s 
current  solution.  A  negative  reduced  cost  indicates  that  xf  should  be  translated  into  a  col¬ 
umn  for  the  RMP,  and  inserted  into  it,  and  a  non-negative  reduced  cost  indicates  that  no 
favorable  column  currently  exists  for  facility  i.  (A  positive  reduced  cost  for  subproblem 
i  can  arise  if  K,  contains  only  the  null  schedule.) 
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3.3.  Solving  the  column-generation  subproblems 

The  subproblems  CSUB;(jt,  fij )  are  multi-dimensional  knapsack  problems  (Weingart- 
ner  and  Ness  [58])  with  elastic  penalties  in  each  dimension;  Kleywegt  et  al.  [34]  refer 
to  these  as  static  stochastic  knapsack  problems.  We  solve  them  through  straightforward 
branch  and  bound,  except  that  we  add  “explicit  constraint  branching”  (Appleget  and 
Wood  [3])  by  defining  the  general  integer  variables  gj  and  adding  the  following  con¬ 
straint  to  each  subproblem  i: 


y ~]xjj  -  gj  =  0.  (33) 

jeJ 

The  variable  gj  is  an  “ECB  variable”  and  receives  a  higher  priority  for  branching  than 
does  any  Xjj .  Intuitively,  constraint  branching  provides  a  better  balanced  branch-and- 
bound  enumeration  tree,  and  this  tends  to  reduce  enumeration  (see  Ryan  and  Foster 
[48]). 

3.4.  Solving  the  LP-relaxation  of  the  master  problem 

Branch-and-price  algorithms  (e.g.,  Savelsbergh  [49],  Barnhart  et  al.  [7],  Silva  [52])  are 
appearing  as  complements  to  the  branch-and-cut  algorithms  that  are  commonly  imple¬ 
mented  in  commercial  MIP  solvers.  B&P  combines  a  branch-and-bound  algorithm  with 
a  column-generation  procedure.  Achieving  good  performance  with  column-generation 
is  difficult  (Liibbecke  and  Desrosiers  [41]),  but  a  number  of  enhancements  to  the  basic 
procedure  can  help.  “Duals  stabilization”  comprises  the  most  important  enhancement  in 
our  experience,  so  we  describe  that  here  briefly.  (See  du  Merle  et  al.  [27]  and  Silva  [52] 
for  more  detail.) 

Duals  stabilization  attempts  to  accelerate  the  column-generation  process  that  solves 
CSFLP’s  LP  relaxation.  We  follow  du  Merle  et  al.  [27]  for  this  purpose,  and  incorporate 
an  elastic  dynamic  trust  region  for  dual  variables.  The  trust  region  is  always  centered  on 
the  most  recent  dual  solution.  It  is  elastic  because  penalized  violation  of  the  nominal  trust 
region  is  allowed,  and  it  is  dynamic  because  its  width  and  penalties  are  adjusted  contin¬ 
ually.  This  trust-region  mechanism  is  implemented  by  turning  master-problem  equality 
constraints  into  elastic  ranged  constraints.  The  primal  (master-problem)  elastic  penalties 
define  the  dual  trust  region’s  limits,  while  primal  ranges  define  the  dual  penalties,  i.e., 
the  penalties  applied  if  the  dual  variables  fall  outside  the  nominal  trust  region. 

A  trust  region  of  some  sort  makes  sense  in  this  context  because  (a)  the  column- 
generation  mechanism,  when  viewed  in  the  dual,  is  essentially  Benders  decomposition 
(Benders  [9]),  and  (b)  the  Benders  master  problem  appears  to  benefit  from  the  use  of 
trust  regions  (e.g..  Brown  et  al.  [13],  Linderoth  and  Wright  [37]).  Of  course,  many 
variants  on  trust  regions  could  be  applied  to  our  problems,  but  this  one  is  simple  and 
has  proven  effective  in  recent  column-generation  experiments  (Silva  [52],  Singh  et  al. 
[53]).  (Actually,  experiments  in  [53]  suggest  an  even  better  method  for  “duals  stabiliza¬ 
tion,”  one  which  we  have  not  tested:  Simply  use  interior-point  duals  as  provided  by  an 
interior-point  solution  of  LP-RMR) 
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3.5.  Computational  results 


We  implement  B&P  using  software  from  the  Computational  INfrastructure  for  Opera¬ 
tions  Research  (“COIN-OR,”  or  simply  “COIN”),  which  provides  a  repository  of  distinct 
libraries  that  can  be  integrated  to  build  optimization  algorithms  (Lougee-Heimer  [39]). 
The  COIN  library  labeled  “BCP”  provides  the  basic  framework  for  a  B&P  algorithm 
(Ralphs  and  Ladanyi  [45]).  This  library’s  design  anticipates  a  parallel/distributed  envi¬ 
ronment,  and,  unfortunately,  the  protocol  that  emulates  this  environment  in  our  serial 
environment  incurs  some  computational  overhead.  This  overhead  could  be  avoided  with 
some  additional  programming,  so  the  total  solution  times  reported  here,  denoted  (TT), 
exclude  that  overhead.  However,  we  note  that  the  true  CPU  times  for  our  implementa¬ 
tions  never  exceed  TT  by  more  than  10%,  and  the  mean  overhead  for  all  problems  is 
only  3.1%. 

We  have  implemented  our  B&P  algorithm  using  COIN’S  open  solver  interface  (OSI), 
coupled  with  CPLEX  8.0.  The  linear  relaxation  of  the  RMP  and  the  subproblems  are 
submitted  to  CPLEX’s  LP  solver  and  MIP  solver,  respectively.  We  carry  out  these  tests 
on  a  networked  workstation,  a  Dell  Dimension  340  with  a  2  GHz,  Pentium  IV  processor 
and  1  GB  of  RAM.  For  comparison,  we  also  directly  solve  the  extensive  formulations 
of  SFLP  using  CPLEX  8.0,  and  report  these  solution  times  under  “IP”  in  the  Tables  1 
and  2. 

We  note  that  the  binary  status  of  the  variables  xtj  in  the  extensive  formulation  would 
force  continuous  versions  of  the  variables  xj  to  be  binary  in  any  optimal  solution.  Con¬ 
sequently,  we  need  not  indicate  to  the  solver  that  the  xj  are  binary.  However,  we  find  that 
solution  times  are  shorter,  on  average,  when  we  specify  the  xj  to  be  binary  and  set  the 
branching  priorities  on  these  variables  to  be  higher  than  those  on  the  x/j .  Thus,  we  force 
branching  first  on  the  important  decisions,  i.e.,  whether  or  not  to  open  a  plant,  rather 
than  on  the  relatively  less  important,  and  more  numerous,  customer-to-plant  assignment 
decisions. 


We  investigate  eight  groups  of  problems.  Each  group  is  defined  by  problem  size, 
meaning  “number  of  facilities-number  of  customers,”  and  these  sizes  are:  5-15,  5- 
30,  8-24,  8-48,  10-30,  10-40,  10-50  and  10-60.  For  each  problem  size,  we  consider 
instances  with  1,  10  or  50  scenarios.  (We  consider  larger  problems  with  more  scenarios 
later.)  Because  run  times  vary  somewhat  between  randomly  generated  instances  of  the 
same  size,  we  examine  five  different  instances  for  each  combination  of  problem  size  and 
number  of  scenarios.  All  problems  in  this  paper  are  solved  to  optimality. 

To  generate  the  test  problems,  we  first  create  a  reference  problem  -  the  superscript  “R” 
below  stands  for  “reference”  -  according  to  the  following  rules:  (a)  Customer  demands 
df  are  integers  taken  from  a  discrete  uniform  distribution  17(5,25),  (b)  transportation 
costs  cf  are  integers  from  U(  15,25),  (c)  facility  capacities  are  uf  —  0.8 df/\I\, 
and  (d)  the  fixed  costs  are  cj  =  puf  for  some  cost-per-unit-capacity  conversion  con¬ 
stant  p,  which  is  1.5  for  these  examples.  Chu  and  Beasley  [19]  use  rules  (a)  and  (b) 
to  generate  certain  instances  of  the  generalized-assignment  problem.  For  our  stochastic 
instances,  demands  dj  are  uniformly  distributed  integers  within  ±20%  of  df,  capaci¬ 
ties  uj  are  ±10%  of  uf,  and  fixed  costs  are  simply  cj  =  cjR.  Also,  facility  i  pays  an 
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Table  1.  Total  time  (TT)  in  CPU  seconds,  to  solve  randomly  generated  SFLPs  with  scenario  uncertainty.  Three 
different  algorithms  solve  five  problem  instances  for  each  combination  of  size  and  number  of  scenarios.  (Note: 
All  generated  scenario  data  for  the  problem  instance  in  row  r,  for  r  =  1, . . .  ,  10,  have  been  reused  in  row  r  + 
5.  This  accounts  for  apparent  correlations  in  runtimes  as  exemplified  by  rows  5,  10  and  15.)  A  time  given  in 
the  bold  font  indicates  the  fastest  among  the  three  alternative  solution  methods 


Number  Problem  Size  (facilities-customers) 
of 

Scenarios 


5-15 

5-30 

8-24 

8-48 

IP 

B&P 
w/o  Stz 

B&P 
w/  Stz 

IP 

B&P 
w/o  Stz 

B&P 
w/  Stz 

IP 

B&P 
w/o  Stz 

B&P 
w/  Stz 

IP 

B&P 
w/o  Stz 

B&P 
w/  Stz 

1 

3.5 

0.5 

0.8 

2.5 

1.7 

1.7 

2274.9 

i.i 

14.5 

5.8 

5.0 

2.7 

1 

0.1 

0.3 

0.8 

1.5 

4.1 

3.9 

* 

2.0 

2.7 

3.8 

3.4 

2.4 

1 

3.2 

0.6 

2.0 

0.9 

1.4 

1.9 

274.1 

1.6 

1.8 

4.5 

6.4 

2.7 

1 

0.2 

0.3 

0.6 

1.6 

2.1 

1.9 

299.8 

1.3 

4.4 

3.8 

3.2 

2.6 

1 

3.6 

0.5 

1.4 

8.6 

2.4 

2.8 

1.8 

0.8 

0.9 

5076.1 

61.2 

44.6 

10 

1.3 

1.1 

1.6 

4.2 

3.7 

3.7 

881.0 

4.8 

9.1 

142.4 

7.9 

4.2 

10 

1.4 

1.0 

1.4 

5.1 

16.6 

30.5 

2987.5 

4.3 

7.0 

39.7 

6.7 

5.7 

10 

1.0 

0.9 

1.1 

1.2 

2.8 

3.0 

231.3 

4.8 

5.7 

20.6 

8.4 

5.2 

10 

0.4 

0.6 

1.3 

2.4 

3.2 

3.6 

16.4 

2.7 

2.9 

19.8 

6.7 

5.9 

10 

1.0 

0.7 

1.4 

9.4 

3.6 

6.8 

3.9 

1.2 

1.3 

* 

29.7 

25.8 

50 

2.5 

2.3 

3.0 

4.9 

10.1 

10.2 

1637.2 

45.5 

29.6 

455.6 

40.0 

58.5 

50 

3.1 

2.3 

4.4 

11.0 

11.1 

11.1 

5579.0 

8.7 

7.6 

45.7 

20.2 

14.0 

50 

2.8 

2.3 

2.7 

3.4 

7.8 

8.5 

1081.6 

8.0 

7.6 

53.4 

25.0 

15.9 

50 

1.3 

1.4 

3.1 

4.3 

8.6 

10.1 

57.4 

5.8 

6.0 

129.1 

24.3 

20.9 

50 

3.9 

2.7 

3.3 

12.2 

12.1 

10.7 

6.6 

4.2 

4.2 

990.2 

80.1 

114.6 

Legend: 

IP: 

B&P  w/o  Stz: 
B&P  w  Stz: 


CPLEX  MIP  solver  solving  the  extensive-form  SFLP 
Branch-and-price  algorithm  without  duals  stabilization  solving  CSFLP 
Branch-and-price  algorithm  with  duals  stabilization  solving  CSFLP 
Problem  not  solved  to  optimality  in  7,200  CPU  seconds. 


outside  purchase.  Clearly,  these  parameter  settings  are  somewhat  arbitrary,  but  testing 
suggests  that,  over  a  wide  range  of  settings,  B&P  remains  a  superior  method  for  solving 
these  problems,  compared  to  solving  them  directly.  As  evidence,  the  reader  will  see  in 
Tables  1  and  2  that  B&P  can  even  solve  single-scenario  problems  faster  than  they  can 
be  solved  directly.  That  is,  B&P  can  be  faster  even  when  parameter  variances  are  0.  Fur¬ 
thermore,  experiments  with  |S|  >  1  not  reported  here  indicate  that  B&P’s  superiority 
only  increases  as  parameter  variances  increase  beyond  the  values  used  in  this  paper. 

Tables  1  and  2  show  TT  for  each  problem  instance  solved  (a)  in  its  extensive  form  by 
IP,  (b)  as  CSFLP  by  basic  B&P  without  duals  stabilization,  and  (c)  as  CSFLP  by  B&P 
with  duals  stabilization.  Parameter  settings  with  duals  stabilization  are  fixed  for  all  prob¬ 
lems  tested,  and  we  set  an  arbitrary  limit  of  7,200  seconds  on  total  allowed  computation 
time. 


3.6.  Discussion 


Both  Tables  1  and  2  provide  stark  evidence  that  branch  and  price  can  be  vastly  superior 
to  branch  and  bound  for  solving  certain  SMIPs,  and  even  for  certain  deterministic  MIPs 
as  indicated  by  the  results  for  the  single-scenario  problems.  One  key  to  this  superiority 
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Table  2.  Total  time  (TT)  in  CPU  seconds,  to  solve  randomly  generated  SFLPs  with  scenario  uncertainty.  Three 
different  algorithms  solve  five  problem  instances  for  each  combination  of  size  and  number  of  scenarios.  This 
table  explores  how  computation  times  change  as  the  ratio  of  facilities  to  customers  decreases.  A  time  given  in 
the  bold  font  indicates  the  fastest  among  the  three  alternative  solution  methods 


Number 

of 

Scenarios 

Problem  Size  (facilities-customers) 

10-30 

10-40 

10-50 

10-60 

IP 

B&P 
w/o  Stz 

B&P 
w/  Stz 

IP 

B&P 
w/o  Stz 

B&P 
w/  Stz 

IP 

B&P 
w/o  Stz 

B&P 
w/  Stz 

IP 

B&P 
w/o  Stz 

B&P 
w/  Stz 

i 

* 

16.7 

17.0 

15.1 

3.0 

2.2 

* 

67.5 

143.0 

21.5 

6.8 

3.8 

i 

3.5 

1.1 

1.1 

7.8 

2.1 

2.1 

* 

53.8 

50.8 

22.3 

11.0 

5.3 

i 

44.8 

1.0 

1.5 

14.8 

2.8 

2.1 

1657 

9.7 

7.2 

14.8 

8.3 

4.4 

i 

1.4 

1.3 

1.5 

* 

36.2 

10.2 

* 

116.8 

99.1 

13.5 

11.9 

5.2 

i 

* 

5.3 

19.5 

6.2 

3.2 

1.7 

50 

5.5 

4.1 

23.8 

13.4 

4.3 

10 

* 

10.3 

12.9 

2323.4 

6.0 

4.0 

* 

20.6 

13.9 

37.2 

16.7 

8.0 

10 

5.2 

2.4 

1.8 

860.1 

5.2 

3.8 

* 

310.4 

424.2 

3163.1 

30.2 

22.4 

10 

7.0 

3.3 

2.6 

262.6 

5.6 

4.1 

* 

18.0 

14.1 

140.5 

23.1 

14.5 

10 

5.0 

2.0 

2.1 

* 

53.6 

77.2 

* 

19.7 

11.7 

45.6 

14.2 

10.8 

10 

* 

15.0 

29.4 

613.5 

5.1 

11.2 

* 

20.5 

41.2 

30.7 

19.7 

9.6 

50 

* 

16.2 

15.6 

3347.3 

20.6 

14.2 

* 

99.1 

127.2 

154.6 

37.7 

18.8 

50 

12.6 

6.5 

7.1 

1894.9 

14.5 

11.6 

* 

37.9 

22.9 

* 

60.0 

50.7 

50 

51.7 

7.7 

10.7 

573.8 

14.8 

11.3 

* 

112.3 

171.7 

132.3 

47.9 

39.9 

50 

16.6 

6.7 

6.0 

* 

30.8 

27.3 

* 

33.9 

27.9 

97.1 

34.4 

21.5 

50 

* 

26.0 

115.0 

3559.3 

17.3 

12.3 

* 

72.7 

111.5 

213.3 

39.7 

22.2 

Legend:  Same  as  Table  1 


clearly  lies  in  the  tighter  LP  lower  bound  provided  by  the  column-oriented  formulation 
versus  the  extensive  formulation:  see  Tables  3  and  4. 

One  might  be  concerned  that  B&P  requires  so  much  overhead  that  it  could  not  be 
effective  for  small  problems.  However,  the  problems  in  Table  1  have  only  five  or  eight 
potential  facilities,  and  IP  outperforms  B&P  only  in  problems  with  five  facilities,  and 
then  only  by  a  small  amount.  Moreover,  average  solution  times  for  B&P  are  at  least  an 
order  of  magnitude  faster  than  IP,  and  IP  cannot  even  solve  two  of  the  problems  within 
the  time  limit  of  7,200  CPU  seconds.  Even  for  small  problems,  B&P  is  a  good  choice. 

Table  2,  which  covers  problems  with  10  facilities  and  30  to  60  customers,  clearly 
shows  that  B&P  solution  times  are  more  stable  and  suffer  less  than  IP  when  the  num¬ 
ber  of  scenarios  increases.  Observe  that  (a)  23  problem  instances  out  of  60  could  not  be 
solved  by  IP  within  7,200  CPU  seconds,  (b)  IP  never  outperforms  B&P,  and  (c)  B&P  can 
be  orders  of  magnitude  faster  than  solving  the  original  problem,  even  for  single-scenario 
instances,  i.e.,  for  deterministic  problems. 

Table  5  explores  the  computational  limits  of  our  current  B&P  implementation  by 
covering  a  wider  range  of  problem  sizes  and  number  of  samples  than  do  Tables  1  and  2. 
Camm  et  al.  [16]  solve  a  facility-location  model  for  a  commercial  application  with  17 
potential  facilities  and  123  customer  zones,  so  our  largest  problem  is  roughly  the  same 
size  as  at  least  one  real-world  problem.  We  can  see  here  (and  to  a  degree  in  Tables  1  and 
2)  that  solution  times  tend  to  increase  only  slowly,  perhaps  linearly,  with  the  number 
of  scenarios.  Thus,  the  number  of  scenarios  does  not  seem  to  be  a  strongly  limiting 
factor  with  the  B&P  methodology.  This  bodes  well  for  solving  an  SMIP  through  sam¬ 
pled  approximating  problems,  since  the  probability  of  identifying  the  optimal  solution 
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Table  3.  Integrality  gaps  compared  between  the  extensive  formulation  of  SFLP  (SFLP)  and  the  column-ori¬ 
ented  formulation  (CSFLP).  These  results  correspond  to  the  problems  in  Table  1 


Number  Problem  Size  (facilities-customers) 
of 

Scenarios 


5-15 

5-30 

8-24 

8^18 

SFLP  (%)  CSFLP  (%) 

SFLP (%)  CSFLP (%) 

SFLP  (%) 

CSFLP  (%) 

SFLP  (%) 

CSFLP  (%) 

1 

18.64 

0.00 

18.15 

0.00 

25.29 

0.00 

22.18 

0.00 

1 

28.48 

0.00 

22.33 

0.00 

27.25 

0.05 

23.26 

0.00 

1 

48.20 

0.00 

26.65 

0.00 

27.00 

0.00 

24.67 

0.00 

1 

16.78 

0.04 

18.43 

0.00 

34.65 

0.00 

22.65 

0.00 

1 

22.16 

0.00 

19.26 

0.00 

22.27 

0.00 

22.49 

0.05 

10 

19.21 

0.03 

17.11 

0.00 

25.93 

0.07 

23.28 

0.00 

10 

34.97 

0.00 

23.14 

0.04 

26.21 

0.02 

21.50 

0.00 

10 

37.36 

0.00 

24.38 

0.00 

28.52 

0.00 

23.79 

0.03 

10 

18.62 

0.00 

17.77 

0.00 

31.37 

0.00 

23.29 

0.00 

10 

20.09 

0.00 

19.47 

0.00 

20.93 

0.00 

21.13 

0.00 

50 

19.37 

0.00 

17.05 

0.00 

24.75 

0.07 

23.66 

0.01 

50 

37.02 

0.00 

23.15 

0.00 

25.90 

0.00 

21.39 

0.00 

50 

36.55 

0.00 

24.15 

0.00 

28.68 

0.00 

24.03 

0.02 

50 

18.87 

0.00 

17.86 

0.00 

30.04 

0.00 

23.73 

0.00 

50 

20.33 

0.00 

19.22 

0.00 

21.05 

0.00 

20.57 

0.02 

for  a  discrete  stochastic  program  increases  exponentially  with  the  number  of  sampled 
scenarios  (Kleywegt  et  al.  [34]). 

Table  5  also  shows  that  problem  size,  in  terms  of  facilities  and  customers,  is  a  stron¬ 
ger  limiting  factor  in  solving  SFLP  by  B&P.  We  discuss  potential  reasons  for  this  in  the 
following  sections. 

4.  Solving  a  special  case  of  SFLP  exactly 

Here  we  investigate  a  special  case  of  the  SFLP  in  which  uncertain  parameters  are  con¬ 
tinuously  distributed  random  variables  satisfying  certain  independence  requirements. 
This  model  paradigm  appears  frequently  in  the  literature  (e.g.,  Louveaux  and  Peeters 
[40],  Laporte  et  al.  [36]);  however,  such  models  are  rarely  solved  exactly  as  we  shall 
solve  SFLP.  Even  if  the  reader  believes  that  some  of  our  independence  assumptions  are 
unreasonable  in  a  real-world  facility-location  problem  -  we  require  independence  of 
demands,  for  instance,  and  this  seems  particularly  unlikely  in  the  SFLP  -  it  is  instructive 
to  see  that  exact  solutions  can  be  achieved  for  such  a  model  in  the  column-oriented 
framework.  Perhaps  the  independence  assumptions  will  be  more  appropriate  in  other 
applications. 

Consider  a  random  vector  f  =  vec(c,  d.  f )  whose  elements  represent  shipping  costs, 
customer  demands  and  excess-demand  penalties,  respectively;  facility  capacities  u  will 
be  deterministic  initially.  Using  this  definition  of  f ,  and  the  following  detailed  defini¬ 
tions,  the  column-oriented  formulation  for  SFLP  clearly  fits  the  form  of  CSMIP1. 

Data  [units] 

cl  fixed  cost  for  opening  facility  i  [dollars] 

Cij  per-unit  shipping  cost  from  facility  i  to  customer  j  [dollars/ton] 
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Table  4.  Integrality  gaps  compared  between  the  compact  formulation  of  SFLP  (SFLP)  and  the  column-oriented 
formulation  (CSFLP).  These  results  correspond  to  the  problems  in  Table  2 


Number  Problem  Size  (facilities-customers) 
of 

Scenarios 


10-30 

10-40 

10-50 

10-60 

SFLP  (%) 

CSFLP  (%) 

SFLP  (%) 

CSFLP  (%) 

SFLP  (%) 

CSFLP  (%) 

SFLP  (%)  CSFLP  (%) 

1 

22.35 

0.04 

32.07 

0.00 

40.43 

0.07 

20.72 

0.00 

1 

23.00 

0.00 

26.59 

0.00 

41.53 

0.06 

32.15 

0.00 

1 

21.04 

0.00 

26.16 

0.00 

56.96 

0.02 

25.89 

0.00 

1 

23.92 

0.00 

26.59 

0.00 

52.95 

0.10 

23.19 

0.00 

1 

3S.46 

0.08 

27.41 

0.00 

45.97 

0.00 

22.71 

0.00 

10 

22.69 

0.02 

32.30 

0.00 

37.85 

0.02 

21.04 

0.00 

10 

22.28 

0.00 

27.02 

0.00 

38.12 

0.12 

35.53 

0.01 

10 

20.32 

0.00 

25.77 

0.03 

55.77 

0.00 

27.42 

0.02 

10 

21.80 

0.00 

27.55 

0.05 

48.76 

0.04 

23.96 

0.00 

10 

42.51 

0.09 

27.05 

0.00 

46.41 

0.03 

21.74 

0.00 

50 

22.66 

0.00 

33.43 

0.00 

38.46 

0.03 

20.21 

0.00 

50 

21.65 

0.00 

27.64 

0.00 

36.00 

0.00 

35.39 

0.00 

50 

20.73 

0.00 

25.55 

0.00 

55.82 

0.09 

26.89 

0.00 

50 

21.82 

0.00 

28.18 

0.00 

47.38 

0.00 

24.36 

0.00 

50 

42.30 

0.06 

28.03 

0.00 

44.95 

0.06 

22.23 

0.00 

dj  demand  from  customer  j  [tons] 

Cjj  expected  cost  to  supply  all  of  customer /s  demand  from  facility  i,  disregarding 
any  shortfalls  in  facility  i’s  capacity,  i.e.,  djj  =  E[cijdj]  [dollars] 

Ui  capacity  of  facility  i  [tons] 

fi  per-unit  penalty  for  excess  demand  that  must  be  covered  by  facility  i  [dollars/tons] 
Cj  total  expected  cost  of  the  /cth  joint  assignment  of  customers  to  facility  i  [dollars]: 


0 

7  _  \+_ 

c'  +  c,xf  +  E[fj]  E 

{Zdjxfj-UiJ 

if  xf  =  0 
otherwise 


(34) 


where  h+  =  max{0,  h}.  Note  that,  up  to  this  point,  we  only  require  independence  of 
excess-demand  penalties  and  demands. 


4.1.  Normally  distributed  demands 

For  the  special-case  model,  we  assume  that  demands  are  independent  and  normally 
distributed  with  means  m  j  and  variances  vj,  i.e.,  dj  ~  N(mj,  Vj).  For  simplicity,  we 
also  assume  that  all  means  and  variances  are  integral.  The  reader  will  see  that  our 
techniques  easily  extend  to  means  otjm  /  and  variances  PjVj,  where  a/  and  flj  are 
positive  scale  parameters,  and  m  j  and  vj  represent  integers  running  from  0  to  some  finite 
upper  bounds.  Given  the  efficiency  of  the  dynamic-programming  solution  procedure  that 
requires  m  j  and  Vj ,  the  scale  parameters  can  be  quite  small.  Thus,  a  wide  range  of  actual 
mean-and- variance  combinations  can  be  closely  approximated. 

The  RMP  for  this  model  is  identical  to  the  column-oriented  formulation,  CSFLP, 
presented  in  Section  3.2.  To  solve  this  special  case  exactly,  we  will  exactly  solve  the 
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Table  5.  Total  time  (TT)  in  CPU  seconds,  for  randomly  generated  SFLPs  with  scenario  uncertainty.  This  table 
explores  the  computational  limits  of  our  current  B&P  implementation  with  duals  stabilization 

Number 

of 

Scenarios 

Problem  Size  (facilities-customers) 

10-60 

20-60 

15-80 

20-100 

50 

36.4 

64.6 

47.2 

137.2 

50 

35.2 

55.7 

65.1 

257.2 

50 

26.9 

58.1 

54.6 

108.2 

50 

35.9 

62.1 

44.0 

106.0 

50 

26.7 

53.4 

106.6 

154.9 

too 

48.2 

136.8 

132.7 

229.6 

too 

68.4 

96.5 

108.0 

829.6 

too 

56.2 

109.4 

238.3 

165.7 

too 

73.0 

123.9 

75.6 

150.0 

too 

51.0 

90.4 

76.3 

257.0 

200 

134.6 

245.3 

181.5 

513.0 

200 

131.3 

199.4 

302.2 

1294.6 

200 

103.7 

168.9 

762.1 

260.7 

200 

134.0 

184.5 

157.1 

272.2 

200 

103.8 

174.4 

156.5 

555.0 

300 

154.8 

374.4 

305.4 

507.8 

300 

248.7 

305.8 

391.2 

817.0 

300 

144.6 

250.0 

654.1 

535.4 

300 

164.0 

320.9 

493.2 

493.2 

300 

210.8 

276.0 

274.4 

687.1 

subproblems  corresponding  to  the  formulation  (29)— (32).  Given  Equation  (34),  the  sub¬ 
problem  associated  with  facility  i  is  this  “static  stochastic  knapsack  problem”  (Kley  wegt 
et  al.  [34]): 


Zi  =  min 

jt, 76(0,1}  VjeJ 


J2  ~  ni)Xii  +  E[f‘]E 

j£j 


Evaluating  E 


(  E  djXij  -  Ui 
\jeJ 


is  easy,  because  we  know  that 


E[w(m,  i>)+] 


(36) 


for  any  w(m,  u)  ~  N(m,  v),  where  <!>(•)  denotes  the  cumulative  distribution  function 
of  a  standard  normal  random  variable  (see  [34]). 

To  derive  a  complete  solution,  we  place  an  arbitrary  ordering  on  the  elements  of 
J,  i.e.,  J  —  [1, . . .  ,),...  ,  1 7 1 1,  temporarily  ignore  constraint-violation  penalties,  and 
apply  dynamic  programming  to  evaluate  the  functions  g\  (j,  m,  v)  defined  here: 
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j 

giij,  m,  v)  =  min  V  (c,y  -  7r,-/)*fy-, 

XU,  ...  7 

J'=l 

J 

S.t.  /«  ,'Xiy'  <  m 

j'= 1 


(37) 

(38) 


X!  vi'xij'  =  v  ( .39) 

/=i 

Xij'  e  {0,  1)  for/  =  1,...  ,;,  (40) 

for  7?Z  =  0,  ...  ,  //Zniaxi  and  V  =  0,  -  ,  .  ,  Umax*  where  III  max~y^. jaj  m  i  Umax  — 
jeJ  vj  ■  (hi  Practice,  much  smaller  limits  on  mmax  and  i>max  can  and  should  be  used 
for  efficiency’s  sake.)  The  dynamic  program  is  defined  by: 


Initial  Conditions: 

giij,  m,  v)  —  0  forj  =  0,  m  —  0,  i>  =  0,  and 
giij,  m,  v)  =  +oo  for  j  —  0,  m  ^  0,  r/0. 


Recursion: 

gi(j>m ,  v)  =  min  [gi(j  -  l,m,  v),  c,y  -  tvj  +  gi(j  -l,m-  mj ,  u  -  vy)}  (41) 

for  j  =  1,...  ,  |  / 1 ,  771=  1,...  ,  712  max ,  V  =  1 ,  .  .  .  ,  Umax  • 

This  recursion  is  similar  to  that  for  a  two-dimensional  knapsack  problem,  but  for  a  given 
m,  the  objective  value  gj  does  not  depend  on  the  variance  index  v.  This  will  be  used  in 
a  final  calculation,  however. 

Now,  since  the  joint  assignment  x,  yields  an  aggregate  demand  with  distribution 
N^2jeJ  nijXjj,  Y2 /ej  vj*ij^'  the  optimal  objective  value  for  (35)  will  be 

Zj  =  m€(i . mmax},  |g  (|/| ,  m,  v )  +  E[fi]  E  [w(m  -  ut,  v)+]  j  +  a  -  m.  (42) 

ve{l....  ,vmax} 

For  the  case  in  which  facility  capacities  w/  are  also  independent  and  normally  distrib¬ 
uted,  and  independent  from  the  customer  demands,  the  method  just  described  will  work 
after  making  a  single  modification:  The  expectation  E[w(m  —  »,■ ,  u)+]  in  Equation  (42) 
changes  to  E[vu(m  —  E[n;],  v  +  Var  (m,))+]. 


4.2.  Extensions 

The  methods  described  above  will  work  for  problems  with  different  types  of  independent 

random  demand  parameters,  although  numerical  integration  may  be  required.  Suppose, 

~  H  ■ 

for  instance,  that  each  demand  can  be  defined  by  dj  =  rjh  +  mj  where  Hj  is 

a  modest-sized  integer  parameter,  all  rjh  are  independent  and  identically  distributed 
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(iid)  with  mean  0  and  variance  v,  and  all  m ;  are  positive  integers.  Then,  any  aggregate 
demand  djXij  can  be  described  through  its  integral  mean  m ’  =  /gy  rnjxij ,  and 
its  scaled,  integral  variance  v1  —  v^2jeJ  HjXjj.  Thus,  assuming  deterministic  facility 
capacities  Uj,  the  recursion  (41)  applies  with  appropriate  adjustments  for  the  scaled  vari¬ 


ances,  and  E 


djXi 


V 


can  be  computed  by  numerical  integration  over 


the  distribution  of  J ~2jeJ  ^jxU-  '  bis  will  be  straightforward  because  that  distribution  is 
defined  through  the  mean-shifted  convolution  of  t  =  fljcj  xij  iid  random  variables, 
which  is  completely  defined  by  its  bounded  integral  mean  m',  its  bounded,  scaled,  inte¬ 
gral  variance  vr,  and  the  common  distribution  for  rji,. 


4.3.  Computational  results  for  normally  distributed  demands 

To  build  test  instances  here,  we  select  each  of  the  single-scenario  problems  from  Section 
3,  assume  costs  and  penalties  are  the  deterministic  quantities  specified  by  that  problem, 
and  assume  that  each  demand  in  the  problem  represents  the  expected  value  of  an  inde¬ 
pendent,  normally  distributed  demand.  The  variance  for  each  demand  is  then  generated 
as  a  discrete  uniform  random  variable  on  [1,  Vj],  where  V/  is  the  maximum  value  that 
assures  P(dj  <  0)  <  0.001  (e.g.,  Spoerl  and  Wood  [54]).  Table  6  displays  solution 
times  and  integrality  gaps  for  all  40  problem  instances.  We  use  the  software  suite  of 
Section  3  for  solving  these  problems,  but  the  computer  is  an  IBM  G40,  Pentium  IV 
laptop  computer  with  1  GB  of  RAM,  running  at  3  GHz. 


4.4.  Discussion 


As  in  the  Section  3,  we  see  that  duals  stabilization  is  a  useful  enhancement  to  the  B&P 
algorithm. 

With  some  exceptions,  results  displayed  in  Tables  1,  2,  5  and  6  indicate  that  B&P 
performs  better  on  problems  with  a  smaller  customers-to-facilities  ratio.  Similar  results 
have  been  observed  when  solving  generalized-assignment  problems,  where  the  tasks- 
to-agents  ratio  correlates  positively  with  the  number  of  feasible  solutions  the  problem 
instances  have  (Savelsbergh  [49],  Silva  [52]).  In  turn,  the  number  of  feasible  solutions 
correlates  positively  with  the  number  of  columns  in  the  column-oriented  model,  and 
the  more  columns  a  problem  has,  the  harder  it  must  be  to  solve  using  B&P.  The  glaring 
contradiction  to  this  argument  comes  from  the  10-60  column  in  Table  2,  which  indicates 
that  these  problems  are  easier  than  the  10-50  problems.  Theoretically,  the  10-60  prob¬ 
lems  may  have  more  columns  than  the  10-50  problems,  but  the  effective  number,  i.e., 
the  number  of  “cost-effective  columns”  may  be  smaller.  Note  that  Table  4  shows  that  the 
extensive  10-60  models  have  tighter  LP  relaxations  than  do  their  10-50  counterparts, 
which  implies  the  system  is  more  capacity-bound  in  some  sense  (which  is  just  an  artifact 
of  the  problem  generator).  This  may  imply  that  the  only  cost-effective  columns  are  those 
that  use  most  of  a  facility’s  nominal  capacity,  which  are  relatively  few  in  number. 
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Table  6.  Total  solution  time  (TT,  in  CPU  seconds)  to  solve  CSFLP  with  B&P  when  demands  are  independent, 
normal  random  variables.  Each  time  in  the  bold  font  indicates  the  fastest  solution  time  for  the  given  problem. 
(Pentium  IV,  3  GHz  computer  with  1  GB  of  RAM.) 


Problem  Size 

B&P  w/o  Duals  Stabilization 

B&P  with  Duals  Stabilization 

Int.  Gap  (%) 

Facilities 

Customers 

Soln.  Time  (TT,  sec.)  Cols. 

Nodes 

Soln.  Time  (TT. 

sec.)  Cols. 

Nodes 

5 

15 

0.8 

188 

i 

0.4 

194 

13 

0.00 

1.0 

153 

i 

0.5 

181 

1 

0.00 

0.5 

171 

i 

0.4 

143 

1 

0.00 

0.6 

192 

i 

0.7 

205 

1 

0.00 

0.7 

168 

i 

0.6 

191 

1 

0.00 

5 

30 

11.3 

443 

i 

9.7 

404 

1 

0.00 

9.3 

474 

i 

7.3 

437 

1 

0.00 

9.3 

432 

i 

8.8 

438 

1 

0.01 

9.7 

423 

i 

8.4 

416 

1 

0.00 

18.8 

599 

9 

17.5 

695 

7 

0.04 

8 

24 

1.2 

350 

1 

1.1 

337 

1 

0.00 

1.1 

322 

1 

0.9 

310 

1 

0.00 

3.0 

466 

9 

1.7 

326 

3 

0.05 

1.2 

351 

1 

0.8 

296 

1 

0.00 

1.8 

385 

1 

1.6 

363 

1 

0.00 

8 

48 

9.1 

1001 

1 

6.4 

778 

1 

0.00 

10.2 

1006 

1 

7.0 

765 

1 

0.01 

11.2 

1166 

5 

11.5 

1237 

9 

0.02 

10.4 

989 

3 

7.3 

816 

9 

0.01 

11.5 

959 

3 

7.9 

734 

3 

0.00 

10 

30 

3.3 

512 

5 

2.9 

455 

7 

0.01 

2.8 

501 

1 

2.3 

452 

1 

0.02 

2.7 

524 

1 

2.0 

451 

1 

0.00 

3.3 

526 

1 

2.4 

461 

1 

0.01 

2.5 

497 

3 

5.7 

856 

3 

0.03 

10 

40 

6.5 

823 

1 

5.0 

708 

3 

0.00 

9.3 

860 

1 

6.5 

679 

1 

0.00 

13.9 

793 

3 

11.5 

668 

3 

0.02 

14.0 

893 

13 

19.8 

1326 

5 

0.03 

10.0 

791 

1 

6.7 

662 

3 

0.01 

10 

50 

28.5 

1076 

1 

20.7 

862 

3 

0.03 

37.5 

1180 

5 

22.9 

794 

5 

0.01 

48.9 

1624 

17 

65.3 

2144 

19 

0.07 

32.2 

1303 

5 

34.1 

1423 

5 

0.02 

24.5 

1105 

1 

17.9 

837 

1 

0.00 

10 

60 

73.4 

1357 

1 

56.1 

1221 

1 

0.00 

76.1 

1585 

15 

46.2 

1263 

5 

0.01 

71.7 

1434 

1 

50.1 

1213 

1 

0.00 

90.7 

1581 

5 

66.9 

1243 

5 

0.00 

74.0 

1418 

1 

55.6 

1215 

1 

0.00 

5.  Summary,  conclusions  and  recommendations 

5.1.  Summary 

This  paper  proposes  a  column-oriented  model  for  a  class  of  two-stage  stochastic  mixed- 
integer  programs  (SMIPs),  and  describes  examples  of  well-known  deterministic  optimi¬ 
zation  problems  whose  stochastic  versions  fall  into  this  class.  We  show  how  to  solve  such 
problems  with  a  branch-and-price  algorithm  (B&P),  using  a  stochastic  facility-location 
problem  (SFLP)  as  a  computational  example.  We  solve  one  version  with  scenario  uncer¬ 
tainty  as  well  as  one  with  independent,  normally  distributed  parameters.  We  solve  both 
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Table  7.  Larger  instances  of  SFLP.  Total  solution  time  (TT,  in  CPU  seconds)  to  solve  CSFLP  with  B&P  using 
duals  stabilization,  when  demands  are  independent  normal  random  variables.  An  asterisk  indicates  the  problem 
does  not  solve  in  2,400  seconds.  In  this  case,  “Cols.”  is  the  number  of  columns  generated  in  2,400  seconds. 
(Pentium  IV,  3  GHz  computer  with  1  GB  of  RAM.) 


Problem  Size  (facilities-customers) 

10-60 

20-60 

15-80 

20-100 

25-150 

30-200 

Soln.  Time 

Soln.  Time 

Soln.  Time 

Soln.  Time 

Soln.  Time 

Soln.  Time 

(TT,  sec.) 

Cols.  (TT,  sec.) 

Cols. 

(TT,  sec.) 

Cols.  (TT,  sec.) 

Cols. 

(TT,  sec.) 

Cols.  (TT,  sec.) 

Cols. 

56.1 

1221  9.4 

1039 

110.9 

1947  180.2 

2421 

323.2 

3753  1799.7 

5859 

46.2 

1263  9.3 

1068 

80.5 

1719  253.6 

3144 

451.5 

4051  2386.7 

6776 

50.1 

1213  72.9 

4092 

73.8 

1683  117.8 

2284 

1014.4 

5057  * 

(6644) 

66.9 

1243  11.6 

1080 

91.4 

2023  153.0 

2355 

613.6 

4025  * 

(6569) 

55.6 

1215  29.4 

2525 

65.6 

1585  208.9 

2631 

1406.7 

5358  2337.7 

6359 

versions  exactly,  and  demonstrate  how  the  algorithm’s  performance  can  be  improved  by 
“duals  stabilization.’’  The  open-source  code  libraries  of  the  Computational  INfrastruc- 
ture  for  Operations  Research  (COIN-OR)  provide  the  framework  for  our  B&P  algorithm, 
while  CPLEX  8.0  comprises  the  solver  engine. 


5.2.  Conclusions 

This  research  demonstrates  that  B&P  is  an  attractive  and  accessible  method  to  solve 
certain  SMIPs.  For  SFLP  with  scenario  uncertainty,  B&P  can  be  orders  of  magnitude 
faster  than  solving  the  original  problem  by  branch  and  bound,  and  this  can  be  true  even 
for  deterministic,  i.e.,  single-scenario  problems.  Furthermore,  the  ability  to  solve  exactly 
an  SMIP  with  continuously  distributed  parameters  is  highly  unusual  in  the  stochastic¬ 
programming  literature. 


5.3.  Recommendations  for  further  work 

The  B&P  approach  can  be  used  to  solve,  at  least  approximately,  SMIPs  of  the  class 
described  in  Section  2,  but  with  more  general  probability  distributions.  For  instance, 
the  methods  of  “sample-average  approximations”  (Mak  et  al.  [43],  Kleywegt  et  al.  [34]) 
provide  probabilistic  guarantees  on  solution  quality  and  are  based  on  repeated  solutions 
of  sampled  approximating  problems.  However,  a  sampled  approximating  problem  is 
essentially  identical  to  a  stochastic  program  with  scenario  uncertainty. 

Sampled  subproblems  can  be  used  to  identify  favorable  columns  in  a  “nearly  exact 
algorithm,”  too.  Suppose  that  once  a  subproblem’s  integer  variables  are  fixed,  i.e.,  a  col¬ 
umn  of  the  model  has  been  defined,  the  expected  cost  of  that  column  can  be  estimated 
highly  accurately  through  sampling.  This  certainly  holds  for  SFLP,  where  the  penal¬ 
ties  associated  with,  say,  10,000  sampled  demands  for  some  fixed  customers-to-facility 
assignment  can  be  sampled  and  averaged  in  a  fraction  of  a  second.  For  all  intents  and 
purposes  then,  that  average  will  exactly  equal  the  expected  cost  for  the  column,  and 
the  linear-programming  relaxation  of  a  master  problem  containing  such  columns  would 
yield  exact  dual  solutions.  The  only  theoretical  concern  in  this  procedure  is  that  the  solu¬ 
tion  of  a  sampled  subproblem  might  indicate  that  a  column  is  favorable,  but  extended 
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sampling  would  reveal  that  it  is  not.  If  we  solve  many  sampled  subproblems  and  cannot 
identify  an  improving  column,  then  we  might  become  convinced  that  we  have,  in  fact, 
solved  the  LP-RMP.  However,  a  formal  procedure  will  need  to  be  constructed  to  provide 
a  rigorous  “level  of  conviction.” 

Finally,  we  note  that  the  COIN/BCP  software  is  originally  intended  for  use  in  a 
distributed/parallel  environment.  It  will  be  interesting  to  investigate  how  well  our  pro¬ 
cedures  perform  in  such  an  environment. 
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